Enhanced prediction of bolt support drilling pressure using optimized Gaussian process regression

This study introduces a novel method for predicting drilling pressure in bolt support systems by optimizing Gaussian process time series regression (GPR) using hybrid optimization algorithms. The research initially identified significant variations in prediction outcomes based on different kernel functions and historical points combinations in the GPR algorithm. To address this, we explored 160 distinct schemes combining 10 kernel functions and 16 historical points for numerical analysis. Applying three hybrid optimization algorithms—Genetic Algorithm-GPR (GA-GPR), Particle Swarm Optimization-GPR (PSO-GPR), and Ant Colony Algorithm-GPR (ACA-GPR)—we iteratively optimized these key parameters. The PSO-GPR algorithm emerged as the most effective, achieving an 80% prediction accuracy with a deviation range of 1–2 MPa, acceptable in practical drilling operations. This optimization led to the RQ kernel function with 18 historical points as the optimal combination, yielding an RMSE value of 0.0047246, in contrast to the least effective combination (E kernel function with 6 historical points) producing an RMSE of 0.035704. The final outcome of this study is a robust and efficient prediction system for underground bolt support drilling pressure, verified through practical application. This approach significantly enhances the accuracy and efficiency of support systems in geotechnical engineering, demonstrating the practical applicability of the PSO-GPR model in real-world scenarios.

rock at a certain distance may be similar, or there may be local mutations, with its changes being non-linear.The Gaussian process time series regression algorithm demonstrates good adaptability to the time series data analysis of drilling pressure in the process of bolt support.
Gaussian process time series regression is a machine learning method based on Gaussian random processes, kernel functions, Bayesian theory, and time series data input.Compared to other machine learning methods, Gaussian process time series regression is less challenging to implement, can adaptively obtain hyperparameters, flexibly infer non-parameters, and the output value has probabilistic significance.The study of Gaussian process regression originated in the 1940s and was used in signal filtering technology 5 .In recent years, an increasing number of researchers have begun to apply it to practical engineering applications, especially in the fields of astronomy 6 , civil engineering 7 , battery chemistry 8 , and others.
Scholars have contributed significantly to the field of Gaussian process regression in various ways.Yang et al. proposed a deep Gaussian process regression (DGPR) method for evaluating the State of Health (SOH) of lithium-ion batteries, leveraging both Gaussian processes and deep networks 9 .Hong et al. expanded the applicability of the Iterative Power and Amplitude Correction (IPAC) algorithm to simulate non-stationary and non-Gaussian processes, considering five transformation pairs 10 .Zhang et al. successfully reduced computational complexity by exploiting the Kronecker structure presented in the state-space model of spatiotemporal Gaussian processes, verifying their findings through the application of weather data prediction 11 .In the realm of 3D surface modeling, Zhao et al. proposed a spherical multi-output Gaussian process method 12 .Shadab et al. offered a systematic method based on the black box model and experimental design to establish an alternative model for predicting and estimating transformer top oil temperature parameters 13 .Gao et al. introduced a residual fatigue life prediction method for metal materials based on Gaussian process regression, which is employed to predict the residual fatigue life of metal materials under two-step loading 14 .Meanwhile, Zeng et al. utilized Gaussian process regression to predict building power consumption 15 .Jo et al. proposed a machine learning framework for path loss modeling that is based on multi-dimensional regression of artificial neural networks (ANN), variance analysis via Gaussian processes, and feature selection assisted by Principal Component Analysis (PCA) 16 .Lastly, Rong et al. developed a data-driven nonparametric Bayesian model based on Gaussian processes to describe and predict in real time the uncertainties in ship lateral motion and trajectory 17 .
Recent advancements in Gaussian process regression (GPR) have seen its application extend across various domains of geotechnical engineering.For instance, a study applied GPR to predict geological parameters, such as Rock Quality Designation (RQD) along tunnel routes, demonstrating GPR's utility in underground engineering 18 .Similarly, Gaussian process regression has been explored for predicting the California bearing ratio of HARHA-treated expansive soils, highlighting the versatility of GPR in soil mechanics 19 .Additionally, GPR has been employed to estimate the trend and random components of soil properties at different locations, indicating its potential for spatial analysis in geotechnical contexts 20 .Another significant application includes the use of GPR for effective regional post-earthquake building damage inference, showcasing its proficiency in structural damage assessment 21 .Furthermore, the integration of GPR in Monte Carlo simulations for the stability analysis of slopes with varied soil strength underlines its value in risk assessment and safety evaluation in geotechnical engineering 22 .These studies underscore the growing importance and diverse applicability of GPR in geotechnical engineering, yet also indicate the necessity for further exploration, particularly in the area of bolt support drilling pressure prediction, which this study addresses.
The selection of key parameters, such as the kernel function and historical points, in Gaussian process time series regression remains under-explored in geotechnical engineering applications.Existing research often relies on conventional choices like the Squared Exponential (SE) kernel function, with parameters typically based on heuristic methods rather than systematic analysis.This gap highlights a need for a more systematic approach to determine how these critical parameters influence the model's predictive accuracy and reliability in the context of bolt support drilling pressure.
In this study, we propose a refined method for predicting bolt support drilling pressure in underground settings, employing a hybrid optimization algorithm to precisely determine the key parameters in Gaussian process time series regression.Our initial step involves modeling the Gaussian process time series regression algorithm, where we examine the impact of various kernel functions and historical points on prediction outcomes.This is achieved by analyzing 160 different combinations of 10 kernel functions and 16 historical points against underground bolt support drilling pressure data.Subsequently, we utilize genetic, particle swarm, and ant colony optimization algorithms to iteratively refine these parameters.This methodology not only identifies the most effective parameter combinations but also leads to the development of a robust and reliable prediction system for bolt support drilling pressure in underground environments.We further validate the system's practical applicability and generalization capacity by applying the optimized parameters in diverse underground drilling scenarios.
This study introduces a novel approach to predict drilling pressure in bolt support systems, leveraging machine learning techniques.Unlike traditional methods such as the seismic interval velocity method, which predominantly rely on material analysis and have shown limitations in generalization and adaptability, our approach uses a hybrid algorithm to optimize Gaussian Process Time Series Regression parameters.This method not only addresses the gaps in existing literature by providing a more versatile and accurate tool for predicting drilling pressure but also represents the first known application of machine learning for this specific purpose in underground bolt support systems.By combining the strengths of various optimization algorithms (GA-GPR, PSO-GPR, ACA-GPR) and comparing their efficacy, this research offers significant insights into the potential of machine learning in enhancing the predictability and safety of underground drilling operations.

Algorithmic modeling
Building on the established theoretical framework of Gaussian Process Regression (GPR), our approach utilizes this advanced statistical technique to model and analyze drilling pressure data.GPR, a non-parametric, Bayesian approach to regression, offers flexibility and precision, making it ideal for complex datasets like those in geotechnical engineering.The foundational concepts and mathematical formulations of GPR are extensively covered in the seminal work by Rasmussen and Williams 23 .The following Eqs.(1-14) are derived based on their methodology.
Suppose that the sample training set composed of drilling pressure data in the process of roadway excavation is: Consider a sample training set composed of drilling pressure data from the roadway excavation process.Here, x i represents the sample input vector of test data, and y i is the output value.The drilling pressure training input vector matrix and the output values are denoted by X and y vectors respectively.Therefore, the output vector can be expressed as: In this equation, w is the weight vector, assumed to follow a Gaussian distribution w ∼ N(0, p ) , and ε denotes noise, adhering to the standard Gaussian distribution ε ∼ N(0, σ 2 n ) .According to Bayesian theory, the posterior distribution of weight w is represented as: The likelihood p(y | X, w ) is expressed as: The marginal likelihood in formula(3), p(y | X ) , is expressed as: As the marginal likelihood is independent of weight, we can apply the Gaussian distribution's conditional probability formula and incorporate formulae (3), ( 4) and ( 5) to obtain the posterior distribution of weight w .It reveals that the weight w 's posterior distribution follows the Gaussian distribution p(w y, X ) ∼ N(w, A −1 ) , where w = 1 The weight w is then calculated by marginalization, leading to the posterior distribution of the function f * corresponding to the new input sample matrix x * : When handling nonlinear data problems, a mapping function can be employed to map the input vector x into an N dimensional feature space.Here, the training data x corresponds to the input vector Φ(x) , and the training data set X (which includes the training data) corresponds to the input vector matrix �(X) = ϕ (1) N×n .The new input sample x * is associated with the new input sample vector ϕ(x * ) , and the new input sample set X * (which includes the test data) corresponds to the new input sample matrix �(X * ) = ϕ (1) The prediction distribution for nonlinear data problems can then be obtained through: (1) To introduce kernel functions, we propose: where p is a covariance matrix of order N, which is a real symmetric positive semi-definite matrix.Thus, there is an orthogonal matrix U that satisfies p = U U T , where U is an orthogonal matrix composed of p 's eigen- vector columns, and is a diagonal matrix composed of p 's eigenvalues.Hence: Simultaneously solving the above equation with formula (12), we can transform the kernel function (11) into: The prediction distribution of the nonlinear data problem based on the kernel function matrix is: Here, K (X, X) is the kernel function matrix between the input training data, K (X * , X) is the kernel function matrix between the input test data and input training data, and K (X * , X * ) is the kernel function matrix between the input test data.

Algorithm evaluation Indicators
The algorithm's application effect can be quantitatively analyzed by assessing the solution's impact -that is, the predictive effect.Commonly used evaluation indices include the coefficient of determination (R2), mean square error (MSE), root mean square error (RMSE), mean absolute error (MAE), mean relative error (MRE), and the variance of prediction error (PEV), among others.
In this study, due to the large number of samples, we need to consider and evaluate the overall system's solution effect, eliminate the influence of individual maxima or minima on system prediction, and objectively explain the overall prediction performance.This aligns with the evaluation characteristics of RMSE as described in James et al. 24 .
The RMSE evaluation formula is: For the RMSE index, its value range is 0 ~ ∞.The lower the value, the better the evaluation.

Kernel function and hyper-parameter determination
In order to achieve optimum training and prediction results within the Gaussian process for time-series regression, it is vital to select the right kernel function and determine the corresponding hyperparameters.Kernel functions are primarily divided into two categories: isotropic (ISO) kernel functions and automatic relevance determination (ARD) kernel functions.The hyperparameter in the ISO kernel is a one-dimensional scalar, whereas the hyperparameter in the ARD kernel is a multi-dimensional vector, corresponding to the number of input vectors.Therefore, the ARD kernel implies more computational effort, but it also enables the removal of unrelated inputs.In this study, we selected five classical kernel functions from these categories for screening: E, SE, RQ, Matern3/2, Matern5/2 for ISO kernel functions, and ARDE, ARDSE, ARDRQ, ARDMatern3/2, ARDMatern5/2 for ARD kernel functions.
The mathematical formulation and effectiveness of these kernel functions within Gaussian processes have been comprehensively discussed in the work of Schölkopf and Smola.Similarly, the process of hyperparameter determination, especially the computation of the negative log marginal likelihood (Formula 17), is a critical aspect of optimizing the Gaussian process model, as detailed in the same reference 25 .
Taking the Matern3/2 kernel function as an example, its expression involves the sum of unknown parameters σ f and σ l , which are considered hyperparameters.The optimal values are obtained using the method of maximum marginal likelihood.The marginal likelihood function for the input training data set is then the negative logarithm of the marginal likelihood.This function serves as the minimization objective function to search for optimal hyperparameters.

Historical points
Gaussian process time-series regression employs a Gaussian process regression algorithm for time-series data processing.For instance, as time progresses and roadway excavation continues, the drilling pressure value for each supporting bolt can be used to form the input vector for Gaussian process time-series regression.
Assume that during the roadway excavation and support process, the time-series sample vector created by the drilling pressure of the same row of anchor rods is {x 1 , x 2 . . .x n } .The first h group data forms the input vector of the training set, and the h + 1 data acts as the output value of the training set, thus forming the first training sample pair.Then, using the second data to the h + 1 data as the input vector of the training set, and the h + 2 data as the output value of the training set, the second training sample pair is formed, and so on.The number of historical points is crucial to the performance of the algorithm, affecting the number of training sample pairs and the global changes considered by the training samples.

Impact of key parameters on results
We collected drilling pressure test data (Supplementary Information 1) during the excavation process of a roadway, specifically when drilling 1000 mm.The dataset comprised a total of 1000 samples, partitioned into 700 training samples and 300 test samples.We considered 10 kernel functions (E, SE, RQ, Matern3/2, Matern5/2, ARDE, ARDSE, ARDRQ, ARDMatern3/2, ARDMatern5/2) and 16 historical points (5-20), which were combined into 160 distinct schemes for numerical analysis.
We employed the Root Mean Square Error (RMSE) as our evaluation metric.The optimal combination of kernel functions and historical points was found to be the RQ kernel function with 18 historical points, yielding an RMSE value of 0.0047246.The solution curve for this test set is depicted in Fig. 1.Conversely, the least effective combination was the E kernel function with 6 historical points, which produced an RMSE of 0.035704, as shown in Fig. 2.
Comparison of Figs. 1 and 2 reveals considerable variation in calculation results based on the combination of kernel functions and historical points when using the Gaussian process time series algorithm to predict bolt support drilling pressure.The best combination allows for more accurate prediction of drilling pressure trends and a narrower confidence interval, while the worst combination's predicted values diverge significantly from the actual values, offering a wide confidence interval and results with little to no referential value.Thus, it is crucial to select the relatively optimal kernel function and number of historical points before employing the Gaussian process time series algorithm for drilling pressure prediction.( 16) In this section, we used a brute-force approach to compare and select the most effective scheme from the 160 available.While this method provided accurate results, it was time-consuming, requiring approximately 1200 s for the comparative calculation test.Therefore, we investigated the use of an optimization algorithm for iterative optimization of kernel functions and historical points.
Most optimization algorithms, however, necessitate a high number of iterations to achieve accurate optimization goals, making them more suitable for large sample optimization cases.Given our relatively small sample size of kernel function and historical point combinations, and the need to limit the number of iterations due to computation time constraints, identifying an optimization function that aligns with these characteristics is also critical.

Hybrid optimization algorithm screening key parameters
In the process of screening key parameters for our Gaussian Process Regression model, the selection of optimization algorithms plays a pivotal role.We chose Genetic Algorithm (GA), Particle Swarm Optimization (PSO), and Ant Colony Algorithm (ACA) due to their distinct capabilities in navigating complex optimization landscapes.GA is known for its robustness and ability to explore a wide range of solutions, making it ideal for avoiding local optima.PSO is selected for its efficiency and effectiveness in finding optimal solutions quickly, which is crucial for our model's performance.Lastly, ACA is utilized for its unique approach in solving optimization problems, emulating the foraging behavior of ants, which is highly effective in complex scenarios.The amalgamation of these three algorithms allows for a more comprehensive and efficient optimization process, enhancing the accuracy of our model.
In this study, we construct GA-GPR, PSO-GPR, and ACA-GPR hybrid optimization algorithms by combining Genetic Algorithm, Particle Swarm Optimization, Ant Colony Algorithm, and Gaussian Process Regression algorithm.Before operation, a GPR conversion function is constructed with kernel function and historical points as input variables and RMSE as output dependent variables, as shown in Eq. ( 18).
Here, h represents the code of the kernel function, with 1-10 corresponding to 10 types of kernel functions respectively, and i represents the number of historical points (ranging from 5 to 20).Drilling pressure test data (Supplementary Information 1) from a roadway excavation process, drilled to a depth of 1000 mm, is used as the research sample.We use a total of 1000 sample groups, 700 for the training set, and 300 for the test set.The output RMSE gpr is the RMSE evaluation value of the test set corresponding to the Gaussian process time series regression.The ideal end goal of optimization is a test set RMSE with a kernel function = RQ and historical points = 18, i.e., 0.0047246.
All optimization functions aim for the extreme value of Eq. ( 18) as the optimization objective, searching for the optimal combination from 160 combinations of 10 kernel functions and 16 historical points.To highlight the application features of small samples and fewer iterations, all hybrid optimization functions are optimized with a population of 5 and 15 iterations, significantly reducing the total computation amount.While this setting may not converge within a limited number of iterations, it offers a more intuitive comparison of iterative results.

GA-GPR algorithm
The Genetic Algorithm (GA) is an evolutionary computation method that adheres to the principle of "natural selection and survival of the fittest" seen in biology.First proposed by Professor J. Holland of the University of Michigan in 1967, it is an efficient approach for solving optimization problems 26 .In recent years, genetic algorithms have found applications in a variety of fields such as materials 27 , mathematics 28 , and medicine 29 .The genetic algorithm operates in three key steps: (1) Selection, where superior individuals from the current population are chosen.That is, individuals with higher fitness are more likely to contribute offspring to the next generation.(2) Crossover, where a new generation of individuals is created through crossover operations with a certain probability, embodying the concept of information exchange.(3) Mutation, where the value of a string in the string-structured data is randomly altered with a certain probability.
The Genetic Algorithm (GA) operates based on the principles of natural selection and genetic inheritance.Its mechanism of selection, crossover, and mutation effectively explores the solution space, making it particularly suitable for complex optimization problems like ours.GA's ability to avoid local optima and find global solutions enhances the robustness of our GPR model optimization.
In this study, we employed the Genetic Algorithm Toolbox (gatbax) from the University of Sheffield for optimization simulation.We set the selection probability to 0.80, the crossover probability to 0.75, and the mutation probability to 0.02.The selection function includes both roulette wheel selection (rws) and stochastic universal sampling (sus), and there are four types of crossover functions: two-point crossover (xovdp), singlepoint crossover (xovsp), shuffle crossover (xovsh), and multi-point crossover (xovmp).In total, we have eight combination choices for the crossover functions.Each function combination was solved 10 times.Since the gatbax toolbox optimizes for maximum value and the optimization objective of formula ( 18) is minimum value, the optimization fitness function of GA-GPR takes a negative value for formula (18).The final iteration results are presented in Figs. 3,    The evaluation of iterative calculation accuracy suggests that all function combinations reached the target value (− 0.0047246) in fewer than 50% of the 10 GA-GPR simulations, indicating suboptimal performance.When the crossover functions are the same, the accuracy ratio when the selection function equals rws is smaller than when the selection function equals sus.The best solution was found to be the combination of selection function sus and crossover function xovsp, with an accuracy ratio of 40%.
Evaluating the GA-GPR combination with the optimal selection function (sus) and crossover function (xovsp) in terms of iteration calculation time, we found the average calculation time of 10 iterations to be 375 s.This is a reduction of 68.75% compared to the brute-force method, satisfying the requirement for reduced calculation time and demonstrating excellent performance.

PSO-GPR algorithm
Particle Swarm Optimization (PSO) is an optimization algorithm that originated from the study of bird flocking behavior.It was first proposed by Kennedy and Eberhart in 1995 30 .The algorithm is relatively straightforward, updating an individual's position by tracking the individual's best-known position and the best-known position in the swarm 31 , causing all particles to continually gravitate towards the potential optimal solution until convergence 32 .In recent years, the PSO algorithm has been applied to areas such as power systems 33 and task scheduling 34 .
A critical parameter in PSO is the inertia weight (W), which reflects a particle's ability to retain its previous velocity.A larger inertia weight enhances the global search ability, while a smaller one improves local search precision.Therefore, a decreasing inertia weight aligns better with the optimization characteristics.The PSO  algorithm's strength lies in its fast convergence and ability to handle multi-dimensional spaces, contributing significantly to the refinement of our GPR mode.
In this study, we selected four classic decreasing inertia weight schemas and one schema with a constant inertia weight (W = 1), as illustrated in Fig. 11.The optimal fitness function of GA-GPR is formula (19).Moreover, the update speed range for kernel functions and historical points under the PSO-GPR solution is − 4 to 4. The individual speed update proportion parameter (C1) is 1.49445, and the swarm speed update proportion parameter (C2) is 1.49445.The final iteration results are presented in Figs. 12, 13, 14, 15 and 16.Assessing the accuracy of iterative calculations, the ratio of reaching the target value (− 0.0047246) in the 10 PSO-GPR simulations in Figs. 12, 13 and 14 is more than 60%.The ratio for the simulations in Figs. 15 and 16 is 30%.Overall, the performance is satisfactory.The best iteration result is the PSO-GPR simulation when W = 0.9-0.5 * i/15.Under application conditions of a population size of 5 and 15 iterations, the accuracy rate reaches 80%.
In terms of iteration time, the PSO-GPR algorithm (W = 0.9-0.5*i/15)exhibits superior performance.The average calculation time across its 10 iterations is 480 s, which is 60% less than that of the brute-force method.Thus, it satisfactorily meets the requirement of reducing computation time and delivers commendable performance.

ACA-GPR algorithm
The Ant Colony Algorithm (ACA) is an optimization methodology that emulates the foraging behavior of ants in nature 35 .It was first proposed by Marco Dorigo in his doctoral thesis in 1992 36 .As ants forage, they deposit pheromones on their paths, enabling other ants to perceive these pheromones.The concentration of pheromones represents the distance of the path, with a higher concentration indicating a shorter corresponding path 37 .Here, the pheromone can be interpreted as the fitness function value.In recent years, several scholars have proposed enhanced algorithms based on the ACA 38 , mainly applied to scheduling and route planning problems 39 .
Over time, the pheromones left by ants on the path gradually evaporate and disappear.The key parameter here is the pheromone evaporation coefficient (Rho).If Rho is too small, the previously searched path is more likely to be selected again, affecting the global search ability of the algorithm.Conversely, a large Rho can enhance global search ability but may increase many redundant search operations, thereby slowing algorithm convergence.Therefore, Rho's role is akin to inertia weight W in the PSO algorithm, and a decreasing Rho better aligns with optimization characteristics.
Ant Colony Algorithm (ACA) mimics the pheromone trail-laying and following behavior of ants, making it effective in solving discrete optimization problems.It utilizes a probabilistic technique for finding the shortest paths, which is instrumental in exploring the optimal combination of kernel function and historical points in  In this study, we selected four classic Rho decay schemas and one schema with a constant Rho of 0.9, as depicted in Fig. 17.The optimal fitness function of ACA-GPR is also formula (19).Additionally, the transition probability constant of the ACA-GPR solution is 0.1, and the local search step is 2. The final iteration results are presented in Figs. 18, 19, 20, 21 and 22.  Assessing the accuracy of iterative calculations, the ratio of reaching the target value (− 0.0047246) in the 10 ACA-GPR simulations in Figs. 18, 19, 20 and 21 is less than or equal to 30%.The ratio for the simulation in Fig. 22 is 40%.Overall, the performance is subpar.Given limited calculation and iteration times, ACA-GPR exhibits certain shortcomings, such as slow convergence speed and a propensity to get stuck in local optima.
In terms of iteration calculation time, all ACA-GPR simulations exhibit lengthy computation times.Even for the best-performing group (Rho = 0.9), the average time for 10 iterations is 1882s, which is 56.83% longer than the brute-force method.This excess in iteration time undermines the significance of the optimization algorithm, resulting in poor performance.

Result comparison
Taking the drilling pressure test data when drilling 1000 mm in the process of roadway excavation as the research sample, the optimization results of three hybrid optimization algorithms GA-GPR, PSO-GPR, and ACA-GPR on the relative optimal kernel function and the relative optimal historical points are different.Results show that PSO-GPR performs best, followed by GA-GPR, and lastly ACA-GPR.PSO-GPR algorithm is more suitable for screening the key parameter kernel function and historical points in the Gaussian process time series algorithm applied to the prediction of drilling pressure in bolt support.When the sample size is small and the number of iterations is small, PSO-GPR can achieve 80% accuracy with a 60% reduction in time.
In this study, the PSO-GPR model has demonstrated an 80% accuracy rate in predicting the pre-drilling pressure for bolt support.It's important to contextualize this level of accuracy within the operational framework of underground drilling.For practical drilling operations, a deviation of 1-2 MPa in predicted pressure is typically within acceptable limits, considering the use of load-sensitive hydraulic systems in modern drilling equipment.The primary objective of this prediction is not to pinpoint the exact drilling pressure, but rather to enable proactive adjustments, such as reducing drilling speed in advance to protect the drill bit and optimize drilling efficiency.This preemptive adjustment approach, facilitated by our PSO-GPR model, enhances the safety and  In previous studies, the hybrid optimization algorithm is comprehensively evaluated by iterative calculation accuracy and iterative calculation time.Therefore, we introduce the evaluation function of the hybrid optimization algorithm to quantitativelyevaluate the hybrid optimization algorithm.This is shown by Eq. ( 19): In the formula, 'C' signifies the accuracy rate of iterative calculations, and 'T' stands for the reduction rate in computation time compared to the brute-force method.'EF' denotes the evaluation index of the hybrid algorithm, with a higher value indicating superior algorithm performance.

Application of predicting drilling pressure for anchor rod support Methodology
Building on previous research, we have established a predictive system for drilling pressure in underground bolt support.The basic steps are as follows: The application flow chart is shown in Fig. 23.
Step 1: Pressure and displacement sensors are outfitted on the bolt support drill to monitor varying drilling pressures at different depths.Over time, this data is accumulated to form a "drilling depth -drilling pressure" database.
Step 2: The database is partitioned into three sets: the "Optimization Algorithm Determination Set," the "Key Parameter Optimization Set," and the "GPR Learning Set".
Step 3: Using the "Optimization Algorithm Determination Set" as a sample, a quantitative evaluation is conducted using the evaluation function of the hybrid optimization algorithm, allowing the selection of the best optimization algorithm for this data type.The quality of the optimization algorithm is related to the law of data type, which remains constant once determined.
Step 4: Utilizing the "Key Parameter Optimization Set" as a sample, the optimal optimization algorithm is used to select the optimal kernel function and number of historical points for Gaussian Process time series regression of bolt support drilling pressure.This step needs to be performed regularly, updating the optimal parameters.
Step 5: The "GPR Learning Set" serves as the sample for machine learning, utilizing the currently optimal kernel function and number of historical points as parameters.Step 6: Predict the drilling pressure of the subsequent bolt support, and use the predicted pressure value to automatically set the relevant drilling angle, drilling speed, and other parameters.Simultaneously, the actual drilling pressure monitored by the pressure sensor is stored, and the "GPR Learning Set" and the "drilling depth -drilling pressure" database are continuously updated.

Case study
Data for this case study was collected from the feed pressure of the electro-hydraulic control anchor drilling frame used by CCTEG Taiyuan Research Institute Co., Ltd in Shaanxi Huangling No. 2 Coal Mine Co., Ltd.In the rapid excavation system designed and developed by CCTEG Taiyuan Research Institute Co., Ltd., the drill frame can perform one-key automatic drilling, equipped with a displacement sensor and a pressure sensor to monitor and store the corresponding pressure value of the drill frame at different depths in real-time.
The drilling pressure test data at 1000 mm (Supplementary Information 1), 1200 mm (Supplementary Information 2), 2400 mm (Supplementary Information 3), and 3000 mm (Supplementary Information 4) depths during roadway excavation were used for the "Optimization Algorithm Determination Set", the "Key Parameter Optimization Set", and the "GPR Learning Set".The last 300 drilling pressures at 2400 mm and 3000 mm were predicted using the algorithm.
According to the studies, the PSO-GPR (w = 0.9-0.5*i/15)was found to be the best hybrid optimization algorithm for the current data type, as determined through the "Optimization Algorithm Determination Set".Using the "Key Parameter Optimization Set" as input data, the optimal parameters of Gaussian Process time series regression were determined by the PSO-GPR algorithm.The algorithm population number of PSO-GPR (w = 0.9-0.5*i/15)was set to 8, the iteration number to 20, and the remaining parameters were left unchanged.The iterative optimization curve is shown in Fig. 24.
The total time for 20 iterations was 210 s, the optimization result was 0.0316511, the corresponding kernel function was Matern5/2, and the number of historical points was 9.
Using kernel function = Matern5/2 and historical points = 9, the "GPR Learning Set" was learned, and the subsequent drilling pressures were predicted.The test set prediction results for 2400 mm and 3000 mm drilling data are shown in Figs. 25 and 26 respectively.
Figures 25 and 26 reveal excellent prediction results, with a narrow confidence interval bandwidth.The drilling pressure prediction method for bolt support, based on the hybrid optimization algorithm to screen the key parameters of Gaussian process time series regression, exhibits a certain level of generalization.

Methods
In this study, we developed a robust methodology to optimize the Gaussian Process Regression (GPR) model for predicting drilling pressure in underground bolt support, integrating various hybrid algorithms for enhanced accuracy.Central to our optimization approach is the implementation of a fitness function and the application of three distinct optimization algorithms.

GA-GPR algorithm
Utilizing the Genetic Algorithm Toolbox (gatbax) from the University of Sheffield, this algorithm involves genetic operations like selection, crossover, and mutation, targeting the maximization of the fitness value by operating on the negative value of the OptimizeGPR function.The algorithm iteratively evolves solutions, focusing on minimizing the RMSE to identify the optimal historical points and kernel function types.The process is detailed in Fig. 27.

PSO-GPR algorithm
The PSO-GPR algorithm is built upon the principles of particle swarm optimization, where particles' positions and velocities are dynamically adjusted according to their fitness, evaluated by the OptimizeGPR function.The iterative process seeks to efficiently converge to the GPR parameters that minimize RMSE.The algorithm's dynamic adaptation is illustrated in Fig. 28.    1. Gaussian process time series regression, depending on the selection of kernel function and historical points, effectively predicts drilling pressure within a confidence interval.The optimal combination significantly enhances prediction accuracy and narrows the confidence interval.2. Among the evaluated hybrid optimization algorithms (GA-GPR, PSO-GPR, ACA-GPR), PSO-GPR demonstrates superior performance, achieving 80% prediction accuracy with substantial computational efficiency, including a 60% reduction in computation time.This level of accuracy is deemed sufficient for practical applications in underground drilling operations, considering the typical tolerances and operational requirements.3. The established prediction system for underground bolt support proves its generalizability across various test data, highlighting potential improvements in operational efficiency and safety.The 80% accuracy rate of the PSO-GPR model, while not perfect, is adequate for enabling proactive adjustments in drilling operations, such as optimizing drilling speed and angle to enhance equipment safety and drilling efficiency.

Limitations
While the study presents significant advancements, it's important to acknowledge certain limitations.The accuracy of the model, while generally sufficient for operational needs, may not capture all nuances of complex underground environments.Additionally, the model's performance is subject to the quality of available data and may require adjustments when applied to different geological conditions.Future research could focus on enhancing the model's adaptability and testing it in a wider range of scenarios.